# konrad.smolinski@gmail.com
# date: 		26/11/2010
# last update: 	12/05/2011 checked

# info:
#
#
# ------------------------------------------------------------------
# Call for libraries:
rm(list=ls())
source("./rcLibSrc.R")
# ------------------------------------------------------------------
# Set Up
# Parameters
K <- 2 
# set up:
a <- c(0,1)
b <- c(0,1)
d <- c(2,3^(-1/2),3^(-1/2),3^(-1/2))
zval <- c(-0.5,0.5)
R <- length(zval)
xthres <- c(0)
xval <- c(0,1)
K <- length(xval)
R <- length(zval)

#prm <- rcPrm(xval=xval,zval=zval,d=c(1,0.577,0.577),xthres=xthres) # default parameters
#dgp <- rcDgp(prm$a,prm$b,prm$d,prm$xthres,prm$xval,prm$zval)
#print(prm)
# DGP
dgp <- rcDgp(a,b,d,xthres,xval,zval)

# Compute RHS	
rhs <- rcGetRHS(dgp)
rhs

# Identified Set: search through parameter space
# (time consuming; involves integtration)

# GRID
a1seq <- seq(-0.2,0.5,len=6)
a2seq <- seq(0.1,1.6,len=6)
b1seq <- seq(0,0,len=1)
b2seq <- seq(0.1,1.6,len=4)

loIdx <- rcLeaveOutIdx(length(xval))
qs <- smolyak.quad(2,40)

abGrid <- as.matrix(expand.grid(a1seq,a2seq,b1seq,b2seq))
lhs <- rcGetLHS(abGrid,qs,loIdx,xval)
lhs



# ------------------------------------------------------------------
# Save Results:
save(set,dgp,b2vec,a0seq,a1seq,a,b,d,zval,xval,xthres,file='../Data/rcEx2-1-Data.RData')

# plot
polyAlpha <- list()

pdf("../Graphs/rxEx1-2.pdf",width=10,height=6)
par(mfrow=c(2,3))
for(i in c(1,3,4,5,6,7)){
	print(i)
	polyAlpha[[i]] <- ahull(set[[i]],alpha=0.1)
	plot(polyAlpha[[i]],main=paste("b[2]=",b2vec[i],sep="")
	,xlab=expression(a[1])
	,ylab=expression(a[2])
	,wlines='vor',col=c(1,2,0,0,0,2))
	
	abline(v=0,col="red",lty=2)
	abline(h=1,col="red",lty=2)
	points(0,1,pch=19,cex=1.2,col=2)
}
dev.off()